########################################################################
# CREATE MINIMAL DF THAT DESCRIBES TOTAL BARCODES
########################################################################
total_barcodes_plot_df <- r1_scores_df %>%
  group_by(target_site, lab) %>%
  count(target_site, lab,  name = "total_barcodes") %>%
  arrange(target_site, lab)

########################################################################
# PLOT
########################################################################
ggplotly(ggplot(total_barcodes_plot_df) + 
    geom_bar(aes(x = target_site, y = total_barcodes, fill = lab), 
    # dodge prserves vertical position, single groups bars of same metric together
    # stat = identity skips aggregation of values - I provide values
             position = position_dodge(preserve = "single"), stat = "identity") + 
    scale_y_continuous(labels = comma) + # add commas to axis label numbers
    theme_bw() + 
    theme(axis.text.x = element_text(angle = 90)) + 
    labs(x = "Target Site", y = "Total Number of Barcodes", 
         fill = "Lab", main = "Total Number of Barcodes") +
      scale_fill_manual(values=cbPalette))

########################################################################
# TABLE
########################################################################
(total_barcodes_plot_df %>%
    pivot_wider(names_from = "lab",
              values_from = "total_barcodes", 
              values_fill = 0) %>%
  rowwise() %>%
  mutate(AVG = mean(c(NNN03m, NNN04m, NNN05m)),
         SD = sd(c(NNN03m, NNN04m, NNN05m)),
         CV = SD / AVG) %>%
  arrange(desc(CV)) %>%
  kbl() %>%
  kable_styling(bootstrap_options = c("striped", "hover", "condensed", "responsive", html_font = "Cambria")))

########################################################################
# REMOVE ENV VARS 
########################################################################
rm(total_barcodes_plot_df)

####################################################################
## FUNCTION TO MAKE PLOTS
####################################################################
make_qc_reads_plot <- function(df, targetsite){
  plot_df <- df %>%
    filter(target_site == {{targetsite}}) %>% 
    pivot_longer(names_to = "proto_pam", 
                 values_to = "read_count", 
                 cols = c("proto","pam")) %>%
    group_by(target_site)
  
  ## plot 
  ggplotly(ggplot(plot_df) + 
    geom_bar(aes(x = var, y = read_count, fill = lab), 
             position = position_dodge(preserve = "single"), stat = "identity") +
    scale_y_continuous(labels = comma) +
    theme_bw() + 
    theme(axis.text.x = element_text(angle = 90)) + 
    facet_wrap(~proto_pam) + 
    labs(x = "Metric", y = "Read Count", 
         fill = "Lab", title = {{targetsite}}) + scale_fill_manual(values=cbPalette))
}  

####################################################################
## FUNCTION TO MAKE TABLE
####################################################################
make_qc_table <- function(df, targetsite){
  labs <- df$labs
  (plot_df <- df %>%
    filter(target_site == {{targetsite}}) %>% 
    pivot_longer(names_to = "proto_pam", values_to = "read_count", cols = c("proto","pam")) %>%
    group_by(target_site) %>%
    kbl() %>%
    kable_styling(bootstrap_options = c("striped", "hover", "condensed", "responsive", html_font = "Cambria")))
}

make_qc_reads_plot(r1_qc_df, "EMX1")

make_qc_table(r1_qc_df, "EMX1")

make_qc_reads_plot(r1_qc_df, "FANCF")

make_qc_table(r1_qc_df, "FANCF")

make_qc_reads_plot(r1_qc_df, "RNF2")

make_qc_table(r1_qc_df, "RNF2")

############################################################
## CREATE SUMMARY DF 
############################################################
r1_one_seq_score_plot_df <- r1_scores_df %>%
  ## retain only necessary cols, barcode will be identical between two datasets
  select(barcode, oneseq_score, lab, target_site) %>%
  ## list the oneseq score under Karl and NIST-2 accordingly
  pivot_wider(names_from = lab, values_from = oneseq_score, values_fill = NA) %>%
  ## drop any rows that contain NA because those were not found in one of the datasets
  ## comparing only those seqs in common between datasets
  drop_na() %>%
  group_by(target_site)

############################################################
## FUNCTION TO MAKE PLOTS
############################################################
## pass the above df and specify which cols to plot
make_oneseq_scatter_plot <- function(df, lab1, lab2){
  ## reduce df
  plot_df <- df %>%
  select(barcode, target_site, all_of(lab1), all_of(lab2))
  ## PLOT
  plot_df %>%
  ggplot(aes_string(x = {{lab1}}, y = {{lab2}})) + 
  geom_point(alpha = 0.25) + 
  ggpubr::stat_cor(aes(label = ..rr.label..), method = "pearson") + 
  facet_wrap(~target_site, scales = "free_x", ncol = 1) + 
  theme_bw() + 
  labs(x = lab1, 
       y =lab2, 
       title = "ONE-Seq Score of Common Barcodes") 
}

make_oneseq_scatter_plot(r1_one_seq_score_plot_df, "NNN03m", "NNN04m")

make_oneseq_scatter_plot(r1_one_seq_score_plot_df, "NNN04m", "NNN05m")

make_oneseq_scatter_plot(r1_one_seq_score_plot_df, "NNN03m", "NNN05m")

############################################################
## REMOVE OBJECT
############################################################
rm(r1_one_seq_score_plot_df)

## curate summary df to plot normalized pam + proto counts for each lab & target site
## retain only necessary cols, barcode will be identical between two datasets
## drop any rows that contain NA becaue those were not found in one of the datasets
## comparing only those seqs in common between datasets
############################################################
## DF CURATION
############################################################
## PROTO
## EMX1, FANCF, RNF2
normalized_counts_proto_df <- r1_scores_df %>%
  select(barcode, normalized_counts_proto, lab, target_site) %>%
  pivot_wider(names_from = lab, values_from = normalized_counts_proto, values_fill = NA)

## PAM
## EMX1, FACNF, RNF2
normalized_counts_pam_df <- r1_scores_df %>%
  select(barcode, normalized_counts_pam, lab, target_site) %>%
  pivot_wider(names_from = lab, values_from = normalized_counts_pam, values_fill = NA)

############################################################
## FUNCTION
############################################################
########################## PROTO ########################## 
make.proto.norm.plot <- function(df, lab1, lab2){
  plot_df <- df %>%
  select(barcode, target_site, all_of(lab1), all_of(lab2)) %>%
  drop_na() %>%
  group_by(target_site)
  
  ## plot 
  plot_df %>%
  ggplot(aes_string(x = {{lab1}}, y = {{lab2}})) + 
  geom_point(alpha = 0.25) + 
  ggpubr::stat_cor(aes(label = ..rr.label..), method = "pearson") + 
  facet_wrap(~target_site, ncol = 1) + 
  theme_bw() + 
  labs(x = lab1, 
       y = lab2,
       title = "Normalized Counts of Common Barcodes: Proto")
}
########################## PAM ##########################
make.pam.norm.plot <- function(df, lab1, lab2){
  plot_df <- df %>%
  select(barcode, target_site, all_of(lab1), all_of(lab2)) %>%
  drop_na() %>%
  group_by(target_site)
  
  ## plot 
  plot_df %>%
  ggplot(aes_string(x = {{lab1}}, y = {{lab2}})) + 
  geom_point(alpha = 0.25) + 
  ggpubr::stat_cor(aes(label = ..rr.label..), method = "pearson") + 
  facet_wrap(~target_site, ncol = 1) + 
  theme_bw() + 
  labs(x = lab1, 
       y = lab2,
       title = "Normalized Counts of Common Barcodes: PAM")
}

make.proto.norm.plot(normalized_counts_proto_df, "NNN03m", "NNN04m")

make.pam.norm.plot(normalized_counts_pam_df, "NNN03m", "NNN04m")

make.proto.norm.plot(normalized_counts_proto_df, "NNN04m", "NNN05m")

make.pam.norm.plot(normalized_counts_pam_df, "NNN04m", "NNN05m")

make.proto.norm.plot(normalized_counts_proto_df, "NNN03m", "NNN05m")

make.pam.norm.plot(normalized_counts_pam_df, "NNN03m", "NNN05m")

############################################################
## REMOVE OBJECTS
############################################################
rm(normalized_counts_proto_df)
rm(normalized_counts_pam_df)

## curate summary df to plot normalized pam + proto counts for each lab & target site
## retain only necessary cols, barcode will be identical between two datasets
## drop any rows that contain NA becaue those were not found in one of the datasets
## comparing only those seqs in common between datasets
############################################################
## DF CURATION
############################################################

## PROTO
## EMX1, FACNF, RNF2
proto_raw_read_count_df <- r1_scores_df %>%
  select(barcode, proto_raw_read_count, lab, target_site) %>%
  pivot_wider(names_from = lab, values_from = proto_raw_read_count, values_fill = NA)

## PAM
## EMX1, FACNF, RNF2
pam_raw_read_count_df <- r1_scores_df %>%
  select(barcode, pam_raw_read_count, lab, target_site) %>%
  pivot_wider(names_from = lab, values_from = pam_raw_read_count, values_fill = NA)

##########################################################
## FUNCTION
##########################################################
############################## PROTO ##############################
make.proto.raw.counts.plot <- function(df, lab1, lab2){
  
  plot_df <- df %>%
  select(barcode, target_site, lab1, lab2) %>%
  drop_na() %>%
  group_by(target_site)

## plot 
plot_df %>%
  ggplot(aes_string(x = {{lab1}}, y = {{lab2}})) + 
  geom_point(alpha = 0.25) + 
  ggpubr::stat_cor(aes(label = ..rr.label..), method = "pearson") + 
  facet_wrap(~target_site, ncol = 1) + 
  theme_bw() + 
  labs(x = lab1, 
       y = lab2,
       title = "raw Counts of Common Barcodes: Proto")
}
############################## PAM ##############################
make.pam.raw.counts.plot <- function(df, lab1, lab2){
  
  plot_df <- df %>%
  select(barcode, target_site, lab1, lab2) %>%
  drop_na() %>%
  group_by(target_site)

## plot 
plot_df %>%
  ggplot(aes_string(x = {{lab1}}, y = {{lab2}})) + 
  geom_point(alpha = 0.25) + 
  ggpubr::stat_cor(aes(label = ..rr.label..), method = "pearson") + 
  facet_wrap(~target_site, ncol = 1) + 
  theme_bw() + 
  labs(x = lab1, 
       y = lab2,
       title = "raw Counts of Common Barcodes: PAM")
}

make.proto.raw.counts.plot(proto_raw_read_count_df, "NNN03m", "NNN04m")

make.pam.raw.counts.plot(pam_raw_read_count_df, "NNN03m", "NNN04m")

make.proto.raw.counts.plot(proto_raw_read_count_df, "NNN05m", "NNN05m")

make.pam.raw.counts.plot(pam_raw_read_count_df, "NNN05m", "NNN05m")

make.proto.raw.counts.plot(proto_raw_read_count_df, "NNN03m", "NNN05m")

make.pam.raw.counts.plot(pam_raw_read_count_df, "NNN03m", "NNN05m")

##########################################################
## REMOVE OBJECTS
##########################################################
rm(proto_raw_read_count_df)
rm(pam_raw_read_count_df)

########################################################################
# DF CURATION 
########################################################################
r1_venn_barcodes_df <- r1_scores_df %>%
  ## DEFINED AT BEGINNING OF R1 SECTION
  filter(lab %in% labs) %>%
  select(barcode, lab, target_site) %>%
  group_by(target_site, barcode) %>%
  summarise(lab_set = str_c(sort(unique(lab)), collapse = "-"))

########################################################################
# PLOT 
########################################################################
ggplotly(ggplot(data = r1_venn_barcodes_df) + 
  geom_bar(aes(x = fct_rev(as_factor(target_site)), fill = lab_set), 
           # position = "fill"
           position = "fill") + 
  coord_flip() + 
  theme_bw() + 
  labs(x = "Target Site", 
       y = "Proportion of Sites", 
       fill = "Triplicate Combinations", 
       title = "R1 Replicate Barcodes") +
    scale_fill_manual(values=cbPalette))

########################################################################
# TABLE 
########################################################################
(r1_venn_table_df <- r1_venn_barcodes_df  %>% 
  group_by(target_site) %>% 
  count(lab_set)  %>% 
  group_by(target_site) %>% 
  mutate(target_site_total = sum(n)) %>% 
  mutate(pct = round(n/target_site_total * 100, 0),
         count_value = paste0(n, " (", pct, "%)")) %>% 
  select(-n, -pct) %>% 
  pivot_wider(names_from = "lab_set", values_from = "count_value")%>%
  kbl() %>%
  kable_styling(bootstrap_options = c("striped", "hover", "condensed", "responsive", html_font = "Cambria")))

########################################################################
# REMOVE UNNECESSARY OBJECTS FROM ENVIRONMENT 
########################################################################
rm(r1_venn_table_df)

###########################################################################
####### FUNCTION
###########################################################################
  cbPalette <- c("#000000", 
                "#E69F00", 
                "#56B4E9", 
                "#009E73", 
                "#F0E442", 
                "#0072B2", 
                "#D55E00", 
                "#CC79A7",
                "#999999", 
                "#20F03B",
                "#FEC44F", 
                "#D95F0E", 
                "#756BB1", 
                "#FF8A33", 
                "#D7BE07", 
                "#A807D7",
                "dodgerblue2", 
                "#E31A1C", 
                "green4", 
                "#6A3D9A", 
                "#FF7F00", 
                "black", 
                "gold1", 
                "skyblue2",
                "#FB9A99", 
                "palegreen2", 
                "#CAB2D6", 
                "#FDBF6F", 
                "gray70", 
                "khaki2", 
                "maroon", 
                "orchid1",
                "deeppink1",  
                "blue1", 
                "steelblue4", 
                "darkturquoise", 
                "green1", 
                "yellow4", 
                "yellow3",
                "darkorange4", 
                "brown")

######################### PROTO
make.proto.density.plots <- function(df, lablist){
  labs <- c(lablist)
  
  density_plot_labset_df <- df %>%
  filter(lab %in% labs) %>%
  select(barcode, lab, target_site) %>%
  group_by(target_site, barcode) %>%
  summarise(lab_set = str_c(sort(unique(lab)), collapse = "-"))
  
  proto_density_plot_df <- df %>%
  filter(lab %in% labs) %>%
  select(barcode, lab, target_site, proto_raw_read_count) %>%
  right_join(density_plot_labset_df) %>%
  count(target_site, lab_set, proto_raw_read_count, name = "Occurrence")
  
  options(scipen = 999)
  
  plot <- ggplot(proto_density_plot_df) + 
    geom_histogram(aes(x = proto_raw_read_count, fill = lab_set), position = "fill") + 
    scale_x_log10() + 
    facet_wrap(~target_site, scales = "free_y", ncol = 1) +
    theme_bw() +
    scale_fill_brewer(type = "qual") +
    labs(x = "Read Count", y = "Proportion of Sites", fill = "Lab Combinations", main = "Proto Raw Read Counts") +
    scale_fill_manual(values=cbPalette)
  
  return(plot)
}

########################## PAM

make.pam.density.plots <- function(df, lablist){
  labs <- c(lablist)
  
  density_plot_labset_df <- df %>%
  filter(lab %in% labs) %>%
  select(barcode, lab, target_site) %>%
  group_by(target_site, barcode) %>%
  summarise(lab_set = str_c(sort(unique(lab)), collapse = "-"))
  
  pam_density_plot_df <- df %>%
  filter(lab %in% labs) %>%
  select(barcode, lab, target_site, pam_raw_read_count) %>%
  right_join(density_plot_labset_df) %>%
  count(target_site, lab_set, pam_raw_read_count, name = "Occurrence")
  
  options(scipen = 999)
  
  plot <- ggplot(pam_density_plot_df) + 
    geom_histogram(aes(x = pam_raw_read_count, fill = lab_set), position = "fill") + 
    scale_x_log10() + 
    facet_wrap(~target_site, scales = "free_y", ncol = 1) +
    theme_bw() +
    scale_fill_brewer(type = "qual") +
    labs(x = "Read Count", y = "Proportion of Sites", fill = "Lab Combinations", main = "PAM Raw Read Counts") +
    scale_fill_manual(values=cbPalette)
  
  return(plot)
}

lab_list <- c("NNN03m", "NNN04m", "NNN05m")
ggplotly(make.proto.density.plots(r1_scores_df, lab_list))

ggplotly(make.pam.density.plots(r1_scores_df, lab_list))

########################################################################
# CREATE MINIMAL DF THAT DESCRIBES TOTAL BARCODES
########################################################################
total_barcodes_plot_df <- r2_scores_df %>%
  group_by(target_site, lab) %>%
  count(target_site, lab,  name = "total_barcodes") %>%
  arrange(target_site, lab)

########################################################################
# PLOT
######################################################################## 
ggplotly(ggplot(total_barcodes_plot_df) + 
    geom_bar(aes(x = target_site, y = total_barcodes, fill = lab), 
             # dodge prserves vertical position, single groups bars of same metric together
             # stat = identity skips aggregation of values - I provide values
             position = position_dodge(preserve = "single"), stat = "identity") + 
    scale_y_continuous(labels = comma) + # add commas to axis label numbers
    theme_bw() + 
    theme(axis.text.x = element_text(angle = 90)) + 
    labs(x = "Target Site", y = "Total Number of Barcodes", 
         fill = "Lab", main = "Total Number of Barcodes") +
      scale_fill_manual(values=cbPalette))

########################################################################
# TABLE
########################################################################
(total_barcodes_plot_df %>%
    pivot_wider(names_from = "lab",
              values_from = "total_barcodes", 
              values_fill = 0) %>%
  rowwise() %>%
  mutate(AVG = mean(c(NNN06m, NNN07m, NNN08m)),
         SD = sd(c(NNN06m, NNN07m, NNN08m)),
         CV = SD / AVG) %>%
  arrange(desc(CV)) %>%
  kbl() %>%
  kable_styling(bootstrap_options = c("striped", "hover", "condensed", "responsive", html_font = "Cambria")))

########################################################################
# REMOVE ENV VARS 
########################################################################
rm(total_barcodes_plot_df)

## FUNCTIONS ARE IN R1 SECTION
make_qc_reads_plot(r2_qc_df, "EMX1")

make_qc_table(r2_qc_df, "EMX1")

## FUNCTIONS ARE IN R1 SECTION
make_qc_reads_plot(r2_qc_df, "FANCF")

make_qc_table(r2_qc_df, "FANCF")

## FUNCTIONS ARE IN R1 SECTION
make_qc_reads_plot(r2_qc_df, "RNF2")

make_qc_table(r2_qc_df, "RNF2")

############################################################
## CREATE SUMMARY DF 
############################################################
r2_one_seq_score_plot_df <- r2_scores_df %>%
  ## retain only necessary cols, barcode will be identical between two datasets
  select(barcode, oneseq_score, lab, target_site) %>%
  ## list the oneseq score under Karl and NIST-2 accordingly
  pivot_wider(names_from = lab, values_from = oneseq_score, values_fill = NA) %>%
  ## drop any rows that contain NA because those were not found in one of the datasets
  ## comparing only those seqs in common between datasets
  drop_na() %>%
  group_by(target_site)

## FUNCTION IN R1 SECTION
make_oneseq_scatter_plot(r2_one_seq_score_plot_df, "NNN06m", "NNN07m")

## FUNCTION IN R1 SECTION
make_oneseq_scatter_plot(r2_one_seq_score_plot_df, "NNN07m", "NNN08m")

## FUNCTION IN R1 SECTION
make_oneseq_scatter_plot(r2_one_seq_score_plot_df, "NNN06m", "NNN08m")

############################################################
## REMOVE OBJECT
############################################################
rm(r2_one_seq_score_plot_df)

## curate summary df to plot normalized pam + proto counts for each lab & target site
## retain only necessary cols, barcode will be identical between two datasets
## drop any rows that contain NA becaue those were not found in one of the datasets
## comparing only those seqs in common between datasets
############################################################
## DF CURATION
############################################################
################################### PROTO ##################
## EMX1, FANCF, RNF2
normalized_counts_proto_df <- r2_scores_df %>%
  select(barcode, normalized_counts_proto, lab, target_site) %>%
  pivot_wider(names_from = lab, values_from = normalized_counts_proto, values_fill = NA)

################################### PAM ####################
## EMX1, FACNF, RNF2
normalized_counts_pam_df <- r2_scores_df %>%
  select(barcode, normalized_counts_pam, lab, target_site) %>%
  pivot_wider(names_from = lab, values_from = normalized_counts_pam, values_fill = NA)

## FUNCTION DEFINED IN R1 SECTION
make.proto.norm.plot(normalized_counts_proto_df, "NNN06m", "NNN07m")

make.pam.norm.plot(normalized_counts_pam_df, "NNN06m", "NNN07m")

## FUNCTION DEFINED IN R1 SECTION
make.proto.norm.plot(normalized_counts_proto_df, "NNN07m", "NNN08m")

make.pam.norm.plot(normalized_counts_pam_df, "NNN07m", "NNN08m")

## FUNCTION DEFINED IN R1 SECTION
make.proto.norm.plot(normalized_counts_proto_df, "NNN06m", "NNN08m")

make.pam.norm.plot(normalized_counts_pam_df, "NNN06m", "NNN08m")

############################################################
## REMOVE OBJECTS
############################################################
rm(normalized_counts_proto_df)
rm(normalized_counts_pam_df)

## curate summary df to plot normalized pam + proto counts for each lab & target site
## retain only necessary cols, barcode will be identical between two datasets
## drop any rows that contain NA becaue those were not found in one of the datasets
## comparing only those seqs in common between datasets
#############################################
## DF CURATION
#############################################
######################### PROTO
## EMX1, FACNF, RNF2
proto_raw_read_count_df <- r2_scores_df %>%
  select(barcode, proto_raw_read_count, lab, target_site) %>%
  pivot_wider(names_from = lab, values_from = proto_raw_read_count, values_fill = NA)

######################### PAM
## EMX1, FACNF, RNF2
pam_raw_read_count_df <- r2_scores_df %>%
  select(barcode, pam_raw_read_count, lab, target_site) %>%
  pivot_wider(names_from = lab, values_from = pam_raw_read_count, values_fill = NA)

## FUNCTION DEFINED IN R1 SECTION
make.proto.raw.counts.plot(proto_raw_read_count_df, "NNN06m", "NNN07m")

make.pam.raw.counts.plot(pam_raw_read_count_df, "NNN06m", "NNN07m")

## FUNCTION DEFINED IN R1 SECTION
make.proto.raw.counts.plot(proto_raw_read_count_df, "NNN07m", "NNN08m")

make.pam.raw.counts.plot(pam_raw_read_count_df, "NNN07m", "NNN08m")

## FUNCTION DEFINED IN R1 SECTION
make.proto.raw.counts.plot(proto_raw_read_count_df, "NNN06m", "NNN08m")

make.pam.raw.counts.plot(pam_raw_read_count_df, "NNN06m", "NNN08m")

##########################################################
## REMOVE OBJECTS
##########################################################
rm(proto_raw_read_count_df)
rm(pam_raw_read_count_df)

##########################################################
## DF CURATION
##########################################################
r2_venn_barcodes_df <- r2_scores_df %>%
  filter(lab %in% labs) %>%
  select(barcode, lab, target_site) %>%
  group_by(target_site, barcode) %>%
  summarise(lab_set = str_c(sort(unique(lab)), collapse = "-"))

##########################################################
## PLOT
##########################################################
ggplotly(ggplot(data = r2_venn_barcodes_df) + 
  geom_bar(aes(x = fct_rev(as_factor(target_site)), fill = lab_set), 
           # position = "fill"
           position = "fill") + 
  coord_flip() + 
  theme_bw() + 
  labs(x = "Target Site", 
       y = "Proportion of Sites", 
       fill = "Triplicate Combinations", 
       title = "R2 Replicate Barcodes") +
    scale_fill_manual(values=cbPalette))

##########################################################
## TABLE
##########################################################
(r2_venn_table_df <- r2_venn_barcodes_df  %>% 
  group_by(target_site) %>% 
  count(lab_set)  %>% 
  group_by(target_site) %>% 
  mutate(target_site_total = sum(n)) %>% 
  mutate(pct = round(n/target_site_total * 100, 0),
         count_value = paste0(n, " (", pct, "%)")) %>% 
  select(-n, -pct) %>% 
  pivot_wider(names_from = "lab_set", values_from = "count_value")%>%
  kbl() %>%
  kable_styling(bootstrap_options = c("striped", "hover", "condensed", "responsive", html_font = "Cambria")))

########################################################################
# REMOVE UNNECESSARY OBJECTS FROM ENVIRONMENT 
########################################################################
rm(r2_venn_table_df)

## FUNCTION DEFINED IN R1 SECTION
lab_list <- c("NNN06m", "NNN07m", "NNN08m")
ggplotly(make.proto.density.plots(r2_scores_df, lab_list))

## FUNCTION DEFINED IN R1 SECTION
ggplotly(make.pam.density.plots(r2_scores_df, lab_list))

######################################################################
## curate summary df to plot number of off-target sits per target site & lab
r1_total_off_target_sites_plot_df <- r1_combined_df %>%
  group_by(target_site) %>%
  count(target_site, name = "r1_total_off_target_sites") %>%
  arrange(target_site)

r2_total_off_target_sites_plot_df <- r2_combined_df %>%
  group_by(target_site) %>%
  count(target_site, name = "r2_total_off_target_sites") %>%
  arrange(target_site)

r1_r2_total_off_target_sites_df <- r1_total_off_target_sites_plot_df %>%
  left_join(r2_total_off_target_sites_plot_df) %>%
  pivot_longer(cols = c("r1_total_off_target_sites",
                        "r2_total_off_target_sites"), 
                        names_to = "replicate", values_to = "value")%>%
  mutate(replicate = str_remove(replicate, "_total_off_target_sites")) %>%
  mutate(replicate = str_replace(replicate, "r1", "R1"),
         replicate = str_replace(replicate, "r2", "R2"))

##########
## plot ##
##########
ggplotly(ggplot(r1_r2_total_off_target_sites_df) + 
    geom_bar(aes(x = target_site, y = value, fill = replicate), 
             # dodge prserves vertical position, single groups bars of same metric together
             # stat = identity skips aggregation of values - I provide values
             position = position_dodge(preserve = "single"), stat = "identity") + 
    scale_y_continuous(labels = comma) + # add commas to axis label numbers
    theme_bw() + 
    theme(axis.text.x = element_text(angle = 90)) + 
    labs(x = "Target Site", y = "Total Number of Sites", 
         fill = "Lab", main = "Total Number of Barcodes") +
      scale_fill_manual(values=cbPalette))

(r1_r2_total_off_target_sites_df %>%
    rename(total_sites = value) %>%
    pivot_wider(names_from = "replicate",
              values_from = "total_sites", 
              values_fill = 0) %>%
  rowwise() %>%
  mutate(AVG = mean(c(R1, R2)),
         SD = sd(c(R1, R2)),
         CV = SD / AVG) %>%
  arrange(desc(CV)) %>%
  kbl() %>%
  kable_styling(bootstrap_options = c("striped", "hover", "condensed", "responsive", html_font = "Cambria")))

########################################################
## REMOVE UNNECESSARY ENV VARS ##
########################################################
rm(r1_total_off_target_sites_plot_df)
rm(r2_total_off_target_sites_plot_df)

########################################################
## CURATE R1 DF
########################################################
r1_proto_qc_df <- r1_qc_df %>%
  select(-pam) %>%
  pivot_wider(names_from = lab,
              values_from = proto,
              values_fill = NA) %>%
  ## make all 0's NA so it isn't considered in the pmin function
  mutate_if(is.numeric, ~replace(., . == 0, NA)) %>%
  mutate(R1_proto_value = pmin(NNN03m, NNN04m, NNN05m, na.rm=TRUE)) %>%
  select(-NNN03m, -NNN04m, -NNN05m) %>%
  mutate(lab = "R1") %>%
  rename(proto = R1_proto_value)

r1_pam_qc_df <- r1_qc_df %>%
  select(-proto) %>%
  pivot_wider(names_from = lab,
              values_from = pam,
              values_fill = NA) %>%
  ## make all 0's NA so it isn't considered in the pmin function
  mutate_if(is.numeric, ~replace(., . == 0, NA)) %>%
  mutate(R1_pam_value = pmin(NNN03m, NNN04m, NNN05m, na.rm=TRUE)) %>%
  select(-NNN03m, -NNN04m, -NNN05m)%>%
  mutate(lab = "R1") %>%
  rename(pam = R1_pam_value)

r1_combined_qc_df <- r1_proto_qc_df %>%
  left_join(r1_pam_qc_df)

########################################################
## CURATE R2 DF
########################################################
r2_proto_qc_df <- r2_qc_df %>%
  select(-pam) %>%
  pivot_wider(names_from = lab,
              values_from = proto,
              values_fill = NA) %>%
  ## make all 0's NA so it isn't considered in the pmin function
  mutate_if(is.numeric, ~replace(., . == 0, NA)) %>%
  mutate(R2_proto_value = pmin(NNN06m, NNN07m, NNN08m, na.rm=TRUE)) %>%
  select(-NNN06m, -NNN07m, -NNN08m) %>%
  mutate(lab = "R2") %>%
  rename(proto = R2_proto_value)

r2_pam_qc_df <- r2_qc_df %>%
  select(-proto) %>%
  pivot_wider(names_from = lab,
              values_from = pam,
              values_fill = NA) %>%
  ## make all 0's NA so it isn't considered in the pmin function
  mutate_if(is.numeric, ~replace(., . == 0, NA)) %>%
  mutate(R2_pam_value = pmin(NNN06m, NNN07m, NNN08m, na.rm=TRUE)) %>%
  select(-NNN06m, -NNN07m, -NNN08m) %>%
  mutate(lab = "R2") %>%
  rename(pam = R2_pam_value)

r2_combined_qc_df <- r2_proto_qc_df %>%
  left_join(r2_pam_qc_df)  

########################################################
## JOIN R1 R2 DFS
########################################################
r1_r2_qc_df <- plyr::rbind.fill(r1_combined_qc_df, r2_combined_qc_df)
  
########################################################
## REMOVE UNNECESSARY ENV VARS ##
########################################################
rm(r1_proto_qc_df)
rm(r1_pam_qc_df)
rm(r2_proto_qc_df)
rm(r2_pam_qc_df)
rm(r1_combined_qc_df)
rm(r2_combined_qc_df)

make_qc_reads_plot(r1_r2_qc_df, "EMX1")

make_qc_table(r1_r2_qc_df, "EMX1")